MODULE fefmm_LIP1_Fij

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_build_Fij_LIP1

   REAL(REALK), PARAMETER :: T_VALUE_THR = 0.5D0

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Fij_LIP1(FEdim,NFS,FEgrain,t,Fij)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim, NFS
      REAL(REALK),   INTENT(IN)  :: FEgrain, t
      REAL(REALK),   INTENT(OUT) :: Fij(:,:)

      INTEGER(INTK) :: i,k, u,v

      Fij = zero

      DO i = 1, FEdim 
!      i = (FEdim+1)/2
         v = i + FEdim*((NFS-1)/2)  ! V-point

         DO u = 1, NFS*( FEdim -1 )   ! finite elements

            ! 1st LIP
            Fij(u,i)   = Fij(u,i)   + get_F1_LIP1(t,FEgrain,v,u)
            ! 2nd LIP
            Fij(u+1,i) = Fij(u+1,i) + get_F2_LIP1(t,FEgrain,v,u)

         END DO

      END DO

   END SUBROUTINE fefmm_build_Fij_LIP1

!-------------------------------------------------------------------------------

   FUNCTION get_F1_LIP1(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F1_LIP1

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      F1 = -(exp(-t ** 2 * k ** 2 * lambda ** 2)                        &
           - exp(-t ** 2 * lambda ** 2 * (-0.1D1 + k) ** 2))            &
           / t ** 2 / lambda / 0.2D1

      F2 = ((-k / 0.2D1 + 0.1D1 / 0.2D1) * erf(t * k * lambda) + k      &
           * erf(-lambda * t + t * k * lambda) / 0.2D1                  &
           - erf(-lambda * t + t * k * lambda) / 0.2D1) / t             &
           * sqrt(0.3141592654D1)

      get_F1_LIP1 = F1 + F2

   END FUNCTION get_F1_LIP1

!-------------------------------------------------------------------------------

   FUNCTION get_F2_LIP1(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F2_LIP1

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      F1 = (exp(-t ** 2 * k ** 2 * lambda ** 2)                 &
          - exp(-lambda ** 2 * t ** 2 * (-0.1D1 + k) ** 2))     &
          / t ** 2 / lambda / 0.2D1

      F2 = (k * erf(t * k * lambda) / 0.2D1 - k                  &
           * erf(-lambda * t + t * k * lambda) / 0.2D1)          &
           / t * sqrt(0.3141592654D1)

      get_F2_LIP1 = F1 + F2

   END FUNCTION get_F2_LIP1

!-------------------------------------------------------------------------------

END MODULE fefmm_LIP1_Fij
